!--------------------------------------------------------------------------------------------------!
!   CP2K: A general program to perform molecular dynamics simulations                              !
!   Copyright 2000-2025 CP2K developers group <https://cp2k.org>                                   !
!                                                                                                  !
!   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
!--------------------------------------------------------------------------------------------------!

! **************************************************************************************************
!> \brief Initialize a QM/MM calculation
!> \par History
!>      5.2004 created [fawzi]
!> \author Fawzi Mohamed
! **************************************************************************************************
MODULE qmmm_init
   USE atomic_kind_list_types,          ONLY: atomic_kind_list_type
   USE atomic_kind_types,               ONLY: atomic_kind_type,&
                                              get_atomic_kind,&
                                              set_atomic_kind
   USE cell_methods,                    ONLY: read_cell
   USE cell_types,                      ONLY: cell_type,&
                                              use_perd_xyz
   USE cp_control_types,                ONLY: dft_control_type
   USE cp_eri_mme_interface,            ONLY: cp_eri_mme_init_read_input,&
                                              cp_eri_mme_set_params
   USE cp_log_handling,                 ONLY: cp_get_default_logger,&
                                              cp_logger_type
   USE cp_output_handling,              ONLY: cp_print_key_finished_output,&
                                              cp_print_key_unit_nr
   USE cp_subsys_types,                 ONLY: cp_subsys_get,&
                                              cp_subsys_type
   USE cp_units,                        ONLY: cp_unit_from_cp2k,&
                                              cp_unit_to_cp2k
   USE external_potential_types,        ONLY: fist_potential_type,&
                                              get_potential,&
                                              set_potential
   USE force_field_types,               ONLY: input_info_type
   USE force_fields_input,              ONLY: read_gd_section,&
                                              read_gp_section,&
                                              read_lj_section,&
                                              read_wl_section
   USE input_constants,                 ONLY: &
        RADIUS_QMMM_DEFAULT, calc_always, calc_once, do_eri_mme, do_qmmm_gauss, &
        do_qmmm_image_calcmatrix, do_qmmm_image_iter, do_qmmm_link_gho, do_qmmm_link_imomm, &
        do_qmmm_link_pseudo, do_qmmm_pcharge, do_qmmm_swave
   USE input_section_types,             ONLY: section_vals_get,&
                                              section_vals_get_subs_vals,&
                                              section_vals_type,&
                                              section_vals_val_get
   USE kinds,                           ONLY: default_string_length,&
                                              dp
   USE message_passing,                 ONLY: mp_para_env_type
   USE molecule_kind_types,             ONLY: molecule_kind_type
   USE pair_potential_types,            ONLY: pair_potential_reallocate
   USE particle_list_types,             ONLY: particle_list_type
   USE particle_types,                  ONLY: particle_type
   USE pw_env_types,                    ONLY: pw_env_type
   USE qmmm_elpot,                      ONLY: qmmm_potential_init
   USE qmmm_ff_fist,                    ONLY: qmmm_ff_precond_only_qm
   USE qmmm_gaussian_init,              ONLY: qmmm_gaussian_initialize
   USE qmmm_per_elpot,                  ONLY: qmmm_ewald_potential_init,&
                                              qmmm_per_potential_init
   USE qmmm_types_low,                  ONLY: add_set_type,&
                                              add_shell_type,&
                                              create_add_set_type,&
                                              create_add_shell_type,&
                                              qmmm_env_mm_type,&
                                              qmmm_env_qm_type,&
                                              qmmm_links_type
   USE qs_environment_types,            ONLY: get_qs_env,&
                                              qs_environment_type
   USE shell_potential_types,           ONLY: get_shell,&
                                              shell_kind_type
#include "./base/base_uses.f90"

   IMPLICIT NONE
   PRIVATE

   LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qmmm_init'

   PUBLIC :: assign_mm_charges_and_radius, &
             print_qmmm_charges, &
             print_qmmm_links, &
             print_image_charge_info, &
             qmmm_init_gaussian_type, &
             qmmm_init_potential, &
             qmmm_init_periodic_potential, &
             setup_qmmm_vars_qm, &
             setup_qmmm_vars_mm, &
             setup_qmmm_links, &
             move_or_add_atoms, &
             setup_origin_mm_cell

CONTAINS

! **************************************************************************************************
!> \brief Assigns charges and radius to evaluate the MM electrostatic potential
!> \param subsys the subsys containing the MM charges
!> \param charges ...
!> \param mm_atom_chrg ...
!> \param mm_el_pot_radius ...
!> \param mm_el_pot_radius_corr ...
!> \param mm_atom_index ...
!> \param mm_link_atoms ...
!> \param mm_link_scale_factor ...
!> \param added_shells ...
!> \param shell_model ...
!> \par History
!>      06.2004 created [tlaino]
!> \author Teodoro Laino
! **************************************************************************************************
   SUBROUTINE assign_mm_charges_and_radius(subsys, charges, mm_atom_chrg, mm_el_pot_radius, &
                                           mm_el_pot_radius_corr, mm_atom_index, mm_link_atoms, &
                                           mm_link_scale_factor, added_shells, shell_model)
      TYPE(cp_subsys_type), POINTER                      :: subsys
      REAL(KIND=dp), DIMENSION(:), POINTER               :: charges
      REAL(dp), DIMENSION(:), POINTER                    :: mm_atom_chrg, mm_el_pot_radius, &
                                                            mm_el_pot_radius_corr
      INTEGER, DIMENSION(:), POINTER                     :: mm_atom_index, mm_link_atoms
      REAL(dp), DIMENSION(:), POINTER                    :: mm_link_scale_factor
      TYPE(add_shell_type), OPTIONAL, POINTER            :: added_shells
      LOGICAL                                            :: shell_model

      INTEGER                                            :: I, ilink, IndMM, IndShell, ishell
      LOGICAL                                            :: is_shell
      REAL(dp)                                           :: qcore, qi, qshell, rc, ri
      TYPE(atomic_kind_type), POINTER                    :: my_kind
      TYPE(fist_potential_type), POINTER                 :: my_potential
      TYPE(particle_list_type), POINTER                  :: core_particles, particles, &
                                                            shell_particles
      TYPE(particle_type), DIMENSION(:), POINTER         :: core_set, particle_set, shell_set
      TYPE(shell_kind_type), POINTER                     :: shell_kind

      NULLIFY (particle_set, my_kind, added_shells)
      CALL cp_subsys_get(subsys=subsys, particles=particles, core_particles=core_particles, &
                         shell_particles=shell_particles)
      particle_set => particles%els

      IF (ALL(particle_set(:)%shell_index == 0)) THEN
         shell_model = .FALSE.
         CALL create_add_shell_type(added_shells, ndim=0)
      ELSE
         shell_model = .TRUE.
      END IF

      IF (shell_model) THEN
         shell_set => shell_particles%els
         core_set => core_particles%els
         ishell = SIZE(shell_set)
         CALL create_add_shell_type(added_shells, ndim=ishell)
         added_shells%added_particles => shell_set
         added_shells%added_cores => core_set
      END IF

      DO I = 1, SIZE(mm_atom_index)
         IndMM = mm_atom_index(I)
         my_kind => particle_set(IndMM)%atomic_kind
         CALL get_atomic_kind(atomic_kind=my_kind, fist_potential=my_potential, &
                              shell_active=is_shell, shell=shell_kind)
         CALL get_potential(potential=my_potential, &
                            qeff=qi, &
                            qmmm_radius=ri, &
                            qmmm_corr_radius=rc)
         IF (ASSOCIATED(charges)) qi = charges(IndMM)
         mm_atom_chrg(I) = qi
         mm_el_pot_radius(I) = ri
         mm_el_pot_radius_corr(I) = rc
         IF (is_shell) THEN
            IndShell = particle_set(IndMM)%shell_index
            IF (ASSOCIATED(shell_kind)) THEN
               CALL get_shell(shell=shell_kind, &
                              charge_core=qcore, &
                              charge_shell=qshell)
               mm_atom_chrg(I) = qcore
            END IF
            added_shells%mm_core_index(IndShell) = IndMM
            added_shells%mm_core_chrg(IndShell) = qshell
            added_shells%mm_el_pot_radius(Indshell) = ri*1.0_dp
            added_shells%mm_el_pot_radius_corr(Indshell) = rc*1.0_dp
         END IF
      END DO

      IF (ASSOCIATED(mm_link_atoms)) THEN
         DO ilink = 1, SIZE(mm_link_atoms)
            DO i = 1, SIZE(mm_atom_index)
               IF (mm_atom_index(i) == mm_link_atoms(ilink)) EXIT
            END DO
            IndMM = mm_atom_index(I)
            mm_atom_chrg(i) = mm_atom_chrg(i)*mm_link_scale_factor(ilink)
         END DO
      END IF

   END SUBROUTINE assign_mm_charges_and_radius

! **************************************************************************************************
!> \brief Print info on charges generating the qmmm potential..
!> \param mm_atom_index ...
!> \param mm_atom_chrg ...
!> \param mm_el_pot_radius ...
!> \param mm_el_pot_radius_corr ...
!> \param added_charges ...
!> \param added_shells ...
!> \param qmmm_section ...
!> \param nocompatibility ...
!> \param shell_model ...
!> \par History
!>      01.2005 created [tlaino]
!> \author Teodoro Laino
! **************************************************************************************************
   SUBROUTINE print_qmmm_charges(mm_atom_index, mm_atom_chrg, mm_el_pot_radius, mm_el_pot_radius_corr, &
                                 added_charges, added_shells, qmmm_section, nocompatibility, shell_model)
      INTEGER, DIMENSION(:), POINTER                     :: mm_atom_index
      REAL(dp), DIMENSION(:), POINTER                    :: mm_atom_chrg, mm_el_pot_radius, &
                                                            mm_el_pot_radius_corr
      TYPE(add_set_type), POINTER                        :: added_charges
      TYPE(add_shell_type), POINTER                      :: added_shells
      TYPE(section_vals_type), POINTER                   :: qmmm_section
      LOGICAL, INTENT(IN)                                :: nocompatibility, shell_model

      INTEGER                                            :: I, ind1, ind2, IndMM, iw
      REAL(KIND=dp)                                      :: qi, qtot, rc, ri
      TYPE(cp_logger_type), POINTER                      :: logger

      qtot = 0.0_dp
      logger => cp_get_default_logger()
      iw = cp_print_key_unit_nr(logger, qmmm_section, "PRINT%QMMM_CHARGES", &
                                extension=".log")
      IF (iw > 0) THEN
         WRITE (iw, FMT="(/,T2,A)") REPEAT("-", 79)
         WRITE (iw, FMT='(/5X,A)') "MM    POINT CHARGES GENERATING THE QM/MM ELECTROSTATIC POTENTIAL"
         WRITE (iw, FMT="(/,T2,A)") REPEAT("-", 79)
         DO I = 1, SIZE(mm_atom_index)
            IndMM = mm_atom_index(I)
            qi = mm_atom_chrg(I)
            qtot = qtot + qi
            ri = mm_el_pot_radius(I)
            rc = mm_el_pot_radius_corr(I)
            IF (nocompatibility) THEN
               WRITE (iw, '(5X,A9,T15,I5,T28,A8,T38,F12.6,T60,A8,T69,F12.6)') ' MM ATOM:', IndMM, ' RADIUS:', ri, &
                  ' CHARGE:', qi
            ELSE
               WRITE (iw, '(5X,A9,T15,I5,T28,A8,T38,F12.6,T60,A8,T69,F12.6,/,T56,A12,T69,F12.6)') &
                  ' MM ATOM:', IndMM, ' RADIUS:', ri, ' CHARGE:', qi, 'CORR. RADIUS', rc
            END IF
         END DO
         IF (added_charges%num_mm_atoms /= 0) THEN
            WRITE (iw, FMT="(/,T2,A)") REPEAT("-", 79)
            WRITE (iw, '(/5X,A)') "ADDED POINT CHARGES GENERATING THE QM/MM ELECTROSTATIC POTENTIAL"
            WRITE (iw, FMT="(/,T2,A)") REPEAT("-", 79)
            DO I = 1, SIZE(added_charges%mm_atom_index)
               IndMM = added_charges%mm_atom_index(I)
               qi = added_charges%mm_atom_chrg(I)
               qtot = qtot + qi
               ri = added_charges%mm_el_pot_radius(I)
               ind1 = added_charges%add_env(I)%Index1
               ind2 = added_charges%add_env(I)%Index2
               IF (nocompatibility) THEN
                  WRITE (iw, '(5X,A9,I5,T25,A8,T35,F12.6,T50,A8,T59,F12.6,I5,I5)') 'MM POINT:', IndMM, ' RADIUS:', ri, &
                     ' CHARGE:', qi, ind1, ind2
               ELSE
                  WRITE (iw, '(5X,A9,I5,T25,A8,T35,F12.6,T50,A8,T59,F12.6,I5,I5,/,T56,A12,T69,F12.6)') &
                     'MM POINT:', IndMM, ' RADIUS:', ri, ' CHARGE:', qi, ind1, ind2, 'CORR. RADIUS', rc
               END IF
            END DO

         END IF

         IF (shell_model) THEN
            WRITE (iw, FMT="(/,T2,A)") REPEAT("-", 73)
            WRITE (iw, '(/5X,A)') "ADDED SHELL CHARGES GENERATING THE QM/MM ELECTROSTATIC POTENTIAL"
            WRITE (iw, FMT="(/,T2,A)") REPEAT("-", 73)

            DO I = 1, SIZE(added_shells%mm_core_index)
               IndMM = added_shells%mm_core_index(I)
               qi = added_shells%mm_core_chrg(I)
               qtot = qtot + qi
               ri = added_shells%mm_el_pot_radius(I)
               IF (nocompatibility) THEN
                  WRITE (iw, '(7X,A,I5,A8,F12.6,A8,F12.6,3F12.6)') 'SHELL:', IndMM, ' RADIUS:', ri, &
                     ' CHARGE:', qi, added_shells%added_particles(I)%r
               ELSE
                  WRITE (iw, '(7X,A,I5,A8,F12.6,A8,F12.6,A,F12.6)') 'SHELL:', IndMM, ' RADIUS:', ri, &
                     ' CHARGE:', qi, ' CORR. RADIUS', rc
               END IF

            END DO

         END IF

         WRITE (iw, FMT="(/,T2,A)") REPEAT("-", 79)
         WRITE (iw, '(/,T50,A,T69,F12.6)') ' TOTAL CHARGE:', qtot
         WRITE (iw, FMT="(/,T2,A,/)") REPEAT("-", 79)
      END IF
      CALL cp_print_key_finished_output(iw, logger, qmmm_section, &
                                        "PRINT%QMMM_CHARGES")
   END SUBROUTINE print_qmmm_charges

! **************************************************************************************************
!> \brief Print info on qm/mm links
!> \param qmmm_section ...
!> \param qmmm_links ...
!> \par History
!>      01.2005 created [tlaino]
!> \author Teodoro Laino
! **************************************************************************************************
   SUBROUTINE print_qmmm_links(qmmm_section, qmmm_links)
      TYPE(section_vals_type), POINTER                   :: qmmm_section
      TYPE(qmmm_links_type), POINTER                     :: qmmm_links

      INTEGER                                            :: i, iw, mm_index, qm_index
      REAL(KIND=dp)                                      :: alpha
      TYPE(cp_logger_type), POINTER                      :: logger

      logger => cp_get_default_logger()
      iw = cp_print_key_unit_nr(logger, qmmm_section, "PRINT%qmmm_link_info", extension=".log")
      IF (iw > 0) THEN
         IF (ASSOCIATED(qmmm_links)) THEN
            WRITE (iw, FMT="(/,T2, A)") REPEAT("-", 73)
            WRITE (iw, FMT="(/,T31,A)") " QM/MM LINKS "
            WRITE (iw, FMT="(/,T2, A)") REPEAT("-", 73)
            IF (ASSOCIATED(qmmm_links%imomm)) THEN
               WRITE (iw, FMT="(/,T31,A)") " IMOMM TYPE LINK "
               DO I = 1, SIZE(qmmm_links%imomm)
                  qm_index = qmmm_links%imomm(I)%link%qm_index
                  mm_index = qmmm_links%imomm(I)%link%mm_index
                  alpha = qmmm_links%imomm(I)%link%alpha
                  WRITE (iw, FMT="(T2,A,T20,A9,I8,1X,A9,I8,T62,A6,F12.6)") "TYPE: IMOMM", &
                     "QM INDEX:", qm_index, "MM INDEX:", mm_index, "ALPHA:", alpha
               END DO
            END IF
            IF (ASSOCIATED(qmmm_links%pseudo)) THEN
               WRITE (iw, FMT="(/,T31,A)") " PSEUDO TYPE LINK "
               DO I = 1, SIZE(qmmm_links%pseudo)
                  qm_index = qmmm_links%pseudo(I)%link%qm_index
                  mm_index = qmmm_links%pseudo(I)%link%mm_index
                  WRITE (iw, FMT="(T2,A,T20,A9,I8,1X,A9,I8)") "TYPE: PSEUDO", &
                     "QM INDEX:", qm_index, "MM INDEX:", mm_index
               END DO
            END IF
            WRITE (iw, FMT="(/,T2,A,/)") REPEAT("-", 73)
         ELSE
            WRITE (iw, FMT="(/,T2, A)") REPEAT("-", 73)
            WRITE (iw, FMT="(/,T26,A)") " NO QM/MM LINKS DETECTED"
            WRITE (iw, FMT="(/,T2, A)") REPEAT("-", 73)
         END IF
      END IF
      CALL cp_print_key_finished_output(iw, logger, qmmm_section, &
                                        "PRINT%qmmm_link_info")
   END SUBROUTINE print_qmmm_links

! **************************************************************************************************
!> \brief ...
!> \param qmmm_env_qm ...
!> \param para_env ...
!> \param mm_atom_chrg ...
!> \param qs_env ...
!> \param added_charges ...
!> \param added_shells ...
!> \param print_section ...
!> \param qmmm_section ...
!> \par History
!>      1.2005 created [tlaino]
!> \author Teodoro Laino
! **************************************************************************************************
   SUBROUTINE qmmm_init_gaussian_type(qmmm_env_qm, para_env, &
                                      mm_atom_chrg, qs_env, added_charges, added_shells, &
                                      print_section, qmmm_section)
      TYPE(qmmm_env_qm_type), POINTER                    :: qmmm_env_qm
      TYPE(mp_para_env_type), POINTER                    :: para_env
      REAL(KIND=dp), DIMENSION(:), POINTER               :: mm_atom_chrg
      TYPE(qs_environment_type), POINTER                 :: qs_env
      TYPE(add_set_type), POINTER                        :: added_charges
      TYPE(add_shell_type), POINTER                      :: added_shells
      TYPE(section_vals_type), POINTER                   :: print_section, qmmm_section

      INTEGER                                            :: i
      REAL(KIND=dp)                                      :: maxchrg
      REAL(KIND=dp), DIMENSION(:), POINTER               :: maxradius, maxradius2
      TYPE(pw_env_type), POINTER                         :: pw_env

      NULLIFY (maxradius, maxradius2, pw_env)

      maxchrg = MAXVAL(ABS(mm_atom_chrg(:)))
      CALL get_qs_env(qs_env, pw_env=pw_env)
      IF (qmmm_env_qm%add_mm_charges) maxchrg = MAX(maxchrg, MAXVAL(ABS(added_charges%mm_atom_chrg(:))))
      CALL qmmm_gaussian_initialize(qmmm_gaussian_fns=qmmm_env_qm%pgfs, &
                                    para_env=para_env, &
                                    pw_env=pw_env, &
                                    mm_el_pot_radius=qmmm_env_qm%mm_el_pot_radius, &
                                    mm_el_pot_radius_corr=qmmm_env_qm%mm_el_pot_radius_corr, &
                                    qmmm_coupl_type=qmmm_env_qm%qmmm_coupl_type, &
                                    eps_mm_rspace=qmmm_env_qm%eps_mm_rspace, &
                                    maxradius=maxradius, &
                                    maxchrg=maxchrg, &
                                    compatibility=qmmm_env_qm%compatibility, &
                                    print_section=print_section, &
                                    qmmm_section=qmmm_section)

      IF (qmmm_env_qm%move_mm_charges .OR. qmmm_env_qm%add_mm_charges) THEN
         CALL qmmm_gaussian_initialize(qmmm_gaussian_fns=added_charges%pgfs, &
                                       para_env=para_env, &
                                       pw_env=pw_env, &
                                       mm_el_pot_radius=added_charges%mm_el_pot_radius, &
                                       mm_el_pot_radius_corr=added_charges%mm_el_pot_radius_corr, &
                                       qmmm_coupl_type=qmmm_env_qm%qmmm_coupl_type, &
                                       eps_mm_rspace=qmmm_env_qm%eps_mm_rspace, &
                                       maxradius=maxradius2, &
                                       maxchrg=maxchrg, &
                                       compatibility=qmmm_env_qm%compatibility, &
                                       print_section=print_section, &
                                       qmmm_section=qmmm_section)

         SELECT CASE (qmmm_env_qm%qmmm_coupl_type)
         CASE (do_qmmm_gauss, do_qmmm_swave, do_qmmm_pcharge)
            DO i = 1, SIZE(maxradius)
               maxradius(i) = MAX(maxradius(i), maxradius2(i))
            END DO
         END SELECT

         IF (ASSOCIATED(maxradius2)) DEALLOCATE (maxradius2)
      END IF

      IF (qmmm_env_qm%added_shells%num_mm_atoms > 0) THEN

         maxchrg = MAXVAL(ABS(added_shells%mm_core_chrg(:)))
         CALL qmmm_gaussian_initialize(qmmm_gaussian_fns=added_shells%pgfs, &
                                       para_env=para_env, &
                                       pw_env=pw_env, &
                                       mm_el_pot_radius=added_shells%mm_el_pot_radius, &
                                       mm_el_pot_radius_corr=added_shells%mm_el_pot_radius_corr, &
                                       qmmm_coupl_type=qmmm_env_qm%qmmm_coupl_type, &
                                       eps_mm_rspace=qmmm_env_qm%eps_mm_rspace, &
                                       maxradius=maxradius2, &
                                       maxchrg=maxchrg, &
                                       compatibility=qmmm_env_qm%compatibility, &
                                       print_section=print_section, &
                                       qmmm_section=qmmm_section)

         SELECT CASE (qmmm_env_qm%qmmm_coupl_type)
         CASE (do_qmmm_gauss, do_qmmm_swave, do_qmmm_pcharge)
            DO i = 1, SIZE(maxradius)
               maxradius(i) = MAX(maxradius(i), maxradius2(i))
            END DO
         END SELECT

         IF (ASSOCIATED(maxradius2)) DEALLOCATE (maxradius2)

      END IF

      qmmm_env_qm%maxradius => maxradius

   END SUBROUTINE qmmm_init_gaussian_type

! **************************************************************************************************
!> \brief ...
!> \param qmmm_env_qm ...
!> \param mm_cell ...
!> \param added_charges ...
!> \param added_shells ...
!> \param print_section ...
!> \par History
!>      1.2005 created [tlaino]
!> \author Teodoro Laino
! **************************************************************************************************
   SUBROUTINE qmmm_init_potential(qmmm_env_qm, mm_cell, &
                                  added_charges, added_shells, print_section)
      TYPE(qmmm_env_qm_type), POINTER                    :: qmmm_env_qm
      TYPE(cell_type), POINTER                           :: mm_cell
      TYPE(add_set_type), POINTER                        :: added_charges
      TYPE(add_shell_type), POINTER                      :: added_shells
      TYPE(section_vals_type), POINTER                   :: print_section

      CALL qmmm_potential_init(qmmm_coupl_type=qmmm_env_qm%qmmm_coupl_type, &
                               mm_el_pot_radius=qmmm_env_qm%mm_el_pot_radius, &
                               potentials=qmmm_env_qm%potentials, &
                               pgfs=qmmm_env_qm%pgfs, &
                               mm_cell=mm_cell, &
                               compatibility=qmmm_env_qm%compatibility, &
                               print_section=print_section)

      IF (qmmm_env_qm%move_mm_charges .OR. qmmm_env_qm%add_mm_charges) THEN

         CALL qmmm_potential_init(qmmm_coupl_type=qmmm_env_qm%qmmm_coupl_type, &
                                  mm_el_pot_radius=added_charges%mm_el_pot_radius, &
                                  potentials=added_charges%potentials, &
                                  pgfs=added_charges%pgfs, &
                                  mm_cell=mm_cell, &
                                  compatibility=qmmm_env_qm%compatibility, &
                                  print_section=print_section)
      END IF

      IF (qmmm_env_qm%added_shells%num_mm_atoms > 0) THEN

         CALL qmmm_potential_init(qmmm_coupl_type=qmmm_env_qm%qmmm_coupl_type, &
                                  mm_el_pot_radius=added_shells%mm_el_pot_radius, &
                                  potentials=added_shells%potentials, &
                                  pgfs=added_shells%pgfs, &
                                  mm_cell=mm_cell, &
                                  compatibility=qmmm_env_qm%compatibility, &
                                  print_section=print_section)
      END IF

   END SUBROUTINE qmmm_init_potential

! **************************************************************************************************
!> \brief ...
!> \param qmmm_env_qm ...
!> \param qm_cell_small ...
!> \param mm_cell ...
!> \param para_env ...
!> \param qs_env ...
!> \param added_charges ...
!> \param added_shells ...
!> \param qmmm_periodic ...
!> \param print_section ...
!> \param mm_atom_chrg ...
!> \par History
!>      7.2005 created [tlaino]
!> \author Teodoro Laino
! **************************************************************************************************
   SUBROUTINE qmmm_init_periodic_potential(qmmm_env_qm, qm_cell_small, mm_cell, para_env, qs_env, &
                                           added_charges, added_shells, qmmm_periodic, print_section, mm_atom_chrg)
      TYPE(qmmm_env_qm_type), POINTER                    :: qmmm_env_qm
      TYPE(cell_type), POINTER                           :: qm_cell_small, mm_cell
      TYPE(mp_para_env_type), POINTER                    :: para_env
      TYPE(qs_environment_type), POINTER                 :: qs_env
      TYPE(add_set_type), POINTER                        :: added_charges
      TYPE(add_shell_type), POINTER                      :: added_shells
      TYPE(section_vals_type), POINTER                   :: qmmm_periodic, print_section
      REAL(KIND=dp), DIMENSION(:), POINTER               :: mm_atom_chrg

      REAL(KIND=dp)                                      :: maxchrg
      TYPE(dft_control_type), POINTER                    :: dft_control

      IF (qmmm_env_qm%periodic) THEN

         NULLIFY (dft_control)
         CALL get_qs_env(qs_env, dft_control=dft_control)

         IF (dft_control%qs_control%semi_empirical) THEN
            CPABORT("QM/MM periodic calculations not implemented for semi empirical methods")
         ELSE IF (dft_control%qs_control%dftb) THEN
            CALL qmmm_ewald_potential_init(qmmm_env_qm%ewald_env, qmmm_env_qm%ewald_pw, &
                                           qmmm_coupl_type=qmmm_env_qm%qmmm_coupl_type, mm_cell=mm_cell, &
                                           para_env=para_env, qmmm_periodic=qmmm_periodic, print_section=print_section)
         ELSE IF (dft_control%qs_control%xtb) THEN
            CALL qmmm_ewald_potential_init(qmmm_env_qm%ewald_env, qmmm_env_qm%ewald_pw, &
                                           qmmm_coupl_type=qmmm_env_qm%qmmm_coupl_type, mm_cell=mm_cell, &
                                           para_env=para_env, qmmm_periodic=qmmm_periodic, print_section=print_section)
         ELSE

            ! setup for GPW/GPAW
            maxchrg = MAXVAL(ABS(mm_atom_chrg(:)))
            IF (qmmm_env_qm%add_mm_charges) maxchrg = MAX(maxchrg, MAXVAL(ABS(added_charges%mm_atom_chrg(:))))

            CALL qmmm_per_potential_init(qmmm_coupl_type=qmmm_env_qm%qmmm_coupl_type, &
                                         per_potentials=qmmm_env_qm%per_potentials, &
                                         potentials=qmmm_env_qm%potentials, &
                                         pgfs=qmmm_env_qm%pgfs, &
                                         qm_cell_small=qm_cell_small, &
                                         mm_cell=mm_cell, &
                                         compatibility=qmmm_env_qm%compatibility, &
                                         qmmm_periodic=qmmm_periodic, &
                                         print_section=print_section, &
                                         eps_mm_rspace=qmmm_env_qm%eps_mm_rspace, &
                                         maxchrg=maxchrg, &
                                         ncp=qmmm_env_qm%aug_pools(SIZE(qmmm_env_qm%aug_pools))%pool%pw_grid%npts, &
                                         ncpl=qmmm_env_qm%aug_pools(SIZE(qmmm_env_qm%aug_pools))%pool%pw_grid%npts_local)

            IF (qmmm_env_qm%move_mm_charges .OR. qmmm_env_qm%add_mm_charges) THEN

               CALL qmmm_per_potential_init(qmmm_coupl_type=qmmm_env_qm%qmmm_coupl_type, &
                                            per_potentials=added_charges%per_potentials, &
                                            potentials=added_charges%potentials, &
                                            pgfs=added_charges%pgfs, &
                                            qm_cell_small=qm_cell_small, &
                                            mm_cell=mm_cell, &
                                            compatibility=qmmm_env_qm%compatibility, &
                                            qmmm_periodic=qmmm_periodic, &
                                            print_section=print_section, &
                                            eps_mm_rspace=qmmm_env_qm%eps_mm_rspace, &
                                            maxchrg=maxchrg, &
                                            ncp=qmmm_env_qm%aug_pools(SIZE(qmmm_env_qm%aug_pools))%pool%pw_grid%npts, &
                                            ncpl=qmmm_env_qm%aug_pools(SIZE(qmmm_env_qm%aug_pools))%pool%pw_grid%npts_local)
            END IF

            IF (qmmm_env_qm%added_shells%num_mm_atoms > 0) THEN

               CALL qmmm_per_potential_init(qmmm_coupl_type=qmmm_env_qm%qmmm_coupl_type, &
                                            per_potentials=added_shells%per_potentials, &
                                            potentials=added_shells%potentials, &
                                            pgfs=added_shells%pgfs, &
                                            qm_cell_small=qm_cell_small, &
                                            mm_cell=mm_cell, &
                                            compatibility=qmmm_env_qm%compatibility, &
                                            qmmm_periodic=qmmm_periodic, &
                                            print_section=print_section, &
                                            eps_mm_rspace=qmmm_env_qm%eps_mm_rspace, &
                                            maxchrg=maxchrg, &
                                            ncp=qmmm_env_qm%aug_pools(SIZE(qmmm_env_qm%aug_pools))%pool%pw_grid%npts, &
                                            ncpl=qmmm_env_qm%aug_pools(SIZE(qmmm_env_qm%aug_pools))%pool%pw_grid%npts_local)
            END IF

         END IF

      END IF

   END SUBROUTINE qmmm_init_periodic_potential

! **************************************************************************************************
!> \brief ...
!> \param qmmm_section ...
!> \param qmmm_env ...
!> \param subsys_mm ...
!> \param qm_atom_type ...
!> \param qm_atom_index ...
!> \param mm_atom_index ...
!> \param qm_cell_small ...
!> \param qmmm_coupl_type ...
!> \param eps_mm_rspace ...
!> \param qmmm_link ...
!> \param para_env ...
!> \par History
!>      11.2004 created [tlaino]
!> \author Teodoro Laino
! **************************************************************************************************
   SUBROUTINE setup_qmmm_vars_qm(qmmm_section, qmmm_env, subsys_mm, qm_atom_type, &
                                 qm_atom_index, mm_atom_index, qm_cell_small, qmmm_coupl_type, eps_mm_rspace, &
                                 qmmm_link, para_env)
      TYPE(section_vals_type), POINTER                   :: qmmm_section
      TYPE(qmmm_env_qm_type), POINTER                    :: qmmm_env
      TYPE(cp_subsys_type), POINTER                      :: subsys_mm
      CHARACTER(len=default_string_length), &
         DIMENSION(:), POINTER                           :: qm_atom_type
      INTEGER, DIMENSION(:), POINTER                     :: qm_atom_index, mm_atom_index
      TYPE(cell_type), POINTER                           :: qm_cell_small
      INTEGER, INTENT(OUT)                               :: qmmm_coupl_type
      REAL(KIND=dp), INTENT(OUT)                         :: eps_mm_rspace
      LOGICAL, INTENT(OUT)                               :: qmmm_link
      TYPE(mp_para_env_type), POINTER                    :: para_env

      CHARACTER(len=default_string_length)               :: atmname, mm_atom_kind
      INTEGER                                            :: i, icount, ikind, ikindr, my_type, &
                                                            n_rep_val, nkind, size_mm_system
      INTEGER, DIMENSION(:), POINTER                     :: mm_link_atoms
      LOGICAL                                            :: explicit, is_mm, is_qm
      REAL(KIND=dp)                                      :: tmp_radius, tmp_radius_c
      REAL(KIND=dp), DIMENSION(:), POINTER               :: tmp_sph_cut
      TYPE(atomic_kind_list_type), POINTER               :: atomic_kinds
      TYPE(atomic_kind_type), POINTER                    :: atomic_kind
      TYPE(fist_potential_type), POINTER                 :: fist_potential
      TYPE(section_vals_type), POINTER                   :: cell_section, eri_mme_section, &
                                                            image_charge_section, mm_kinds

      NULLIFY (mm_link_atoms, cell_section, tmp_sph_cut)
      NULLIFY (image_charge_section)
      qmmm_link = .FALSE.

      CALL section_vals_get(qmmm_section, explicit=explicit)
      IF (explicit) THEN
         !
         ! QM_CELL
         !
         cell_section => section_vals_get_subs_vals(qmmm_section, "CELL")
         CALL read_cell(qm_cell_small, qm_cell_small, cell_section=cell_section, &
                        check_for_ref=.FALSE., para_env=para_env)
         qm_cell_small%tag = "CELL_QM"
         CALL section_vals_val_get(qmmm_section, "E_COUPL", i_val=qmmm_coupl_type)
         CALL section_vals_val_get(qmmm_section, "EPS_MM_RSPACE", r_val=eps_mm_rspace)
         CALL section_vals_val_get(qmmm_section, "SPHERICAL_CUTOFF", r_vals=tmp_sph_cut)
         CPASSERT(SIZE(tmp_sph_cut) == 2)
         qmmm_env%spherical_cutoff = tmp_sph_cut
         IF (qmmm_env%spherical_cutoff(1) <= 0.0_dp) THEN
            qmmm_env%spherical_cutoff(2) = 0.0_dp
         ELSE
            IF (qmmm_env%spherical_cutoff(2) <= 0.0_dp) qmmm_env%spherical_cutoff(2) = EPSILON(0.0_dp)
            tmp_radius = qmmm_env%spherical_cutoff(1) - 20.0_dp*qmmm_env%spherical_cutoff(2)
            IF (tmp_radius <= 0.0_dp) &
               CALL cp_abort(__LOCATION__, &
                             "SPHERICAL_CUTOFF(1) > 20*SPHERICAL_CUTOFF(1)! Please correct parameters for "// &
                             "the Spherical Cutoff in order to satisfy the previous condition!")
         END IF
         !
         ! Initialization of arrays and core_charge_radius...
         !
         tmp_radius = 0.0_dp
         CALL cp_subsys_get(subsys=subsys_mm, atomic_kinds=atomic_kinds)
         DO Ikind = 1, SIZE(atomic_kinds%els)
            atomic_kind => atomic_kinds%els(Ikind)
            CALL get_atomic_kind(atomic_kind=atomic_kind, &
                                 fist_potential=fist_potential)
            CALL set_potential(potential=fist_potential, &
                               qmmm_radius=tmp_radius, &
                               qmmm_corr_radius=tmp_radius)
            CALL set_atomic_kind(atomic_kind=atomic_kind, &
                                 fist_potential=fist_potential)
         END DO
         CALL setup_qm_atom_list(qmmm_section=qmmm_section, &
                                 qm_atom_index=qm_atom_index, &
                                 qm_atom_type=qm_atom_type, &
                                 mm_link_atoms=mm_link_atoms, &
                                 qmmm_link=qmmm_link)
         !
         ! MM_KINDS
         !
         mm_kinds => section_vals_get_subs_vals(qmmm_section, "MM_KIND")
         CALL section_vals_get(mm_kinds, explicit=explicit, n_repetition=nkind)
         !
         ! Default
         !
         tmp_radius = cp_unit_to_cp2k(RADIUS_QMMM_DEFAULT, "angstrom")
         Set_Radius_Pot_0: DO IkindR = 1, SIZE(atomic_kinds%els)
            atomic_kind => atomic_kinds%els(IkindR)
            CALL get_atomic_kind(atomic_kind=atomic_kind, name=atmname)
            CALL get_atomic_kind(atomic_kind=atomic_kind, &
                                 fist_potential=fist_potential)
            CALL set_potential(potential=fist_potential, qmmm_radius=tmp_radius, &
                               qmmm_corr_radius=tmp_radius)
            CALL set_atomic_kind(atomic_kind=atomic_kind, &
                                 fist_potential=fist_potential)
         END DO Set_Radius_Pot_0
         !
         ! If present overwrite the kind specified in input file...
         !
         IF (explicit) THEN
            DO ikind = 1, nkind
               CALL section_vals_val_get(mm_kinds, "_SECTION_PARAMETERS_", i_rep_section=ikind, &
                                         c_val=mm_atom_kind)
               CALL section_vals_val_get(mm_kinds, "RADIUS", i_rep_section=ikind, r_val=tmp_radius)
               tmp_radius_c = tmp_radius
               CALL section_vals_val_get(mm_kinds, "CORR_RADIUS", i_rep_section=ikind, n_rep_val=n_rep_val)
               IF (n_rep_val == 1) CALL section_vals_val_get(mm_kinds, "CORR_RADIUS", i_rep_section=ikind, &
                                                             r_val=tmp_radius_c)
               Set_Radius_Pot_1: DO IkindR = 1, SIZE(atomic_kinds%els)
                  atomic_kind => atomic_kinds%els(IkindR)
                  CALL get_atomic_kind(atomic_kind=atomic_kind, name=atmname)
                  is_qm = qmmm_ff_precond_only_qm(atmname)
                  IF (TRIM(mm_atom_kind) == atmname) THEN
                     CALL get_atomic_kind(atomic_kind=atomic_kind, &
                                          fist_potential=fist_potential)
                     CALL set_potential(potential=fist_potential, &
                                        qmmm_radius=tmp_radius, &
                                        qmmm_corr_radius=tmp_radius_c)
                     CALL set_atomic_kind(atomic_kind=atomic_kind, &
                                          fist_potential=fist_potential)
                  END IF
               END DO Set_Radius_Pot_1
            END DO
         END IF

         !Image charge section

         image_charge_section => section_vals_get_subs_vals(qmmm_section, "IMAGE_CHARGE")
         CALL section_vals_get(image_charge_section, explicit=qmmm_env%image_charge)

      ELSE
         CPABORT("QMMM section not present in input file!")
      END IF
      !
      ! Build MM atoms list
      !
      size_mm_system = SIZE(subsys_mm%particles%els) - SIZE(qm_atom_index)
      IF (qmmm_link .AND. ASSOCIATED(mm_link_atoms)) size_mm_system = size_mm_system + SIZE(mm_link_atoms)
      ALLOCATE (mm_atom_index(size_mm_system))
      icount = 0

      DO i = 1, SIZE(subsys_mm%particles%els)
         is_mm = .TRUE.
         IF (ANY(qm_atom_index == i)) THEN
            is_mm = .FALSE.
         END IF
         IF (ASSOCIATED(mm_link_atoms)) THEN
            IF (ANY(mm_link_atoms == i) .AND. qmmm_link) is_mm = .TRUE.
         END IF
         IF (is_mm) THEN
            icount = icount + 1
            IF (icount <= size_mm_system) mm_atom_index(icount) = i
         END IF
      END DO
      CPASSERT(icount == size_mm_system)
      IF (ASSOCIATED(mm_link_atoms)) THEN
         DEALLOCATE (mm_link_atoms)
      END IF

      ! Build image charge atom list + set up variables
      !
      IF (qmmm_env%image_charge) THEN
         CALL section_vals_val_get(image_charge_section, "MM_ATOM_LIST", &
                                   explicit=explicit)
         IF (explicit) qmmm_env%image_charge_pot%all_mm = .FALSE.

         IF (qmmm_env%image_charge_pot%all_mm) THEN
            qmmm_env%image_charge_pot%image_mm_list => mm_atom_index
         ELSE
            CALL setup_image_atom_list(image_charge_section, qmmm_env, &
                                       qm_atom_index, subsys_mm)
         END IF

         qmmm_env%image_charge_pot%particles_all => subsys_mm%particles%els

         CALL section_vals_val_get(image_charge_section, "EXT_POTENTIAL", &
                                   r_val=qmmm_env%image_charge_pot%V0)
         CALL section_vals_val_get(image_charge_section, "WIDTH", &
                                   r_val=qmmm_env%image_charge_pot%eta)
         CALL section_vals_val_get(image_charge_section, "DETERM_COEFF", &
                                   i_val=my_type)
         SELECT CASE (my_type)
         CASE (do_qmmm_image_calcmatrix)
            qmmm_env%image_charge_pot%coeff_iterative = .FALSE.
         CASE (do_qmmm_image_iter)
            qmmm_env%image_charge_pot%coeff_iterative = .TRUE.
         END SELECT

         CALL section_vals_val_get(image_charge_section, "RESTART_IMAGE_MATRIX", &
                                   l_val=qmmm_env%image_charge_pot%image_restart)

         CALL section_vals_val_get(image_charge_section, "IMAGE_MATRIX_METHOD", &
                                   i_val=qmmm_env%image_charge_pot%image_matrix_method)

         IF (qmmm_env%image_charge_pot%image_matrix_method == do_eri_mme) THEN
            eri_mme_section => section_vals_get_subs_vals(image_charge_section, "ERI_MME")
            CALL cp_eri_mme_init_read_input(eri_mme_section, qmmm_env%image_charge_pot%eri_mme_param)
            CALL cp_eri_mme_set_params(qmmm_env%image_charge_pot%eri_mme_param, &
                                       hmat=qm_cell_small%hmat, is_ortho=qm_cell_small%orthorhombic, &
                                       zet_min=qmmm_env%image_charge_pot%eta, &
                                       zet_max=qmmm_env%image_charge_pot%eta, &
                                       l_max_zet=0, &
                                       l_max=0, &
                                       para_env=para_env)

         END IF
      END IF

   END SUBROUTINE setup_qmmm_vars_qm

! **************************************************************************************************
!> \brief ...
!> \param qmmm_section ...
!> \param qmmm_env ...
!> \param qm_atom_index ...
!> \param mm_link_atoms ...
!> \param mm_link_scale_factor ...
!> \param fist_scale_charge_link ...
!> \param qmmm_coupl_type ...
!> \param qmmm_link ...
!> \par History
!>      12.2004 created [tlaino]
!> \author Teodoro Laino
! **************************************************************************************************
   SUBROUTINE setup_qmmm_vars_mm(qmmm_section, qmmm_env, qm_atom_index, &
                                 mm_link_atoms, mm_link_scale_factor, &
                                 fist_scale_charge_link, qmmm_coupl_type, &
                                 qmmm_link)
      TYPE(section_vals_type), POINTER                   :: qmmm_section
      TYPE(qmmm_env_mm_type), POINTER                    :: qmmm_env
      INTEGER, DIMENSION(:), POINTER                     :: qm_atom_index, mm_link_atoms
      REAL(KIND=dp), DIMENSION(:), POINTER               :: mm_link_scale_factor, &
                                                            fist_scale_charge_link
      INTEGER, INTENT(OUT)                               :: qmmm_coupl_type
      LOGICAL, INTENT(OUT)                               :: qmmm_link

      LOGICAL                                            :: explicit
      TYPE(section_vals_type), POINTER                   :: qmmm_ff_section

      NULLIFY (qmmm_ff_section)
      qmmm_link = .FALSE.
      CALL section_vals_get(qmmm_section, explicit=explicit)
      IF (explicit) THEN
         CALL section_vals_val_get(qmmm_section, "E_COUPL", i_val=qmmm_coupl_type)
         CALL setup_qm_atom_list(qmmm_section, qm_atom_index=qm_atom_index, qmmm_link=qmmm_link, &
                                 mm_link_atoms=mm_link_atoms, mm_link_scale_factor=mm_link_scale_factor, &
                                 fist_scale_charge_link=fist_scale_charge_link)
         !
         ! Do we want to use a different FF for the non-bonded QM/MM interactions?
         !
         qmmm_ff_section => section_vals_get_subs_vals(qmmm_section, "FORCEFIELD")
         CALL section_vals_get(qmmm_ff_section, explicit=qmmm_env%use_qmmm_ff)
         IF (qmmm_env%use_qmmm_ff) THEN
            CALL section_vals_val_get(qmmm_ff_section, "MULTIPLE_POTENTIAL", &
                                      l_val=qmmm_env%multiple_potential)
            CALL read_qmmm_ff_section(qmmm_ff_section, qmmm_env%inp_info)
         END IF
      END IF
   END SUBROUTINE setup_qmmm_vars_mm

! **************************************************************************************************
!> \brief reads information regarding the forcefield specific for the QM/MM
!>      interactions
!> \param qmmm_ff_section ...
!> \param inp_info ...
!> \par History
!>      12.2004 created [tlaino]
!> \author Teodoro Laino
! **************************************************************************************************
   SUBROUTINE read_qmmm_ff_section(qmmm_ff_section, inp_info)
      TYPE(section_vals_type), POINTER                   :: qmmm_ff_section
      TYPE(input_info_type), POINTER                     :: inp_info

      INTEGER                                            :: n_gd, n_gp, n_lj, n_wl, np
      TYPE(section_vals_type), POINTER                   :: gd_section, gp_section, lj_section, &
                                                            wl_section

!
! NONBONDED
!

      lj_section => section_vals_get_subs_vals(qmmm_ff_section, "NONBONDED%LENNARD-JONES")
      wl_section => section_vals_get_subs_vals(qmmm_ff_section, "NONBONDED%WILLIAMS")
      gd_section => section_vals_get_subs_vals(qmmm_ff_section, "NONBONDED%GOODWIN")
      gp_section => section_vals_get_subs_vals(qmmm_ff_section, "NONBONDED%GENPOT")
      CALL section_vals_get(lj_section, n_repetition=n_lj)
      np = n_lj
      IF (n_lj /= 0) THEN
         CALL pair_potential_reallocate(inp_info%nonbonded, 1, np, lj_charmm=.TRUE.)
         CALL read_lj_section(inp_info%nonbonded, lj_section, start=0)
      END IF
      CALL section_vals_get(wl_section, n_repetition=n_wl)
      np = n_lj + n_wl
      IF (n_wl /= 0) THEN
         CALL pair_potential_reallocate(inp_info%nonbonded, 1, np, williams=.TRUE.)
         CALL read_wl_section(inp_info%nonbonded, wl_section, start=n_lj)
      END IF
      CALL section_vals_get(gd_section, n_repetition=n_gd)
      np = n_lj + n_wl + n_gd
      IF (n_gd /= 0) THEN
         CALL pair_potential_reallocate(inp_info%nonbonded, 1, np, goodwin=.TRUE.)
         CALL read_gd_section(inp_info%nonbonded, gd_section, start=n_lj + n_wl)
      END IF
      CALL section_vals_get(gp_section, n_repetition=n_gp)
      np = n_lj + n_wl + n_gd + n_gp
      IF (n_gp /= 0) THEN
         CALL pair_potential_reallocate(inp_info%nonbonded, 1, np, gp=.TRUE.)
         CALL read_gp_section(inp_info%nonbonded, gp_section, start=n_lj + n_wl + n_gd)
      END IF
      !
      ! NONBONDED14
      !
      lj_section => section_vals_get_subs_vals(qmmm_ff_section, "NONBONDED14%LENNARD-JONES")
      wl_section => section_vals_get_subs_vals(qmmm_ff_section, "NONBONDED14%WILLIAMS")
      gd_section => section_vals_get_subs_vals(qmmm_ff_section, "NONBONDED14%GOODWIN")
      gp_section => section_vals_get_subs_vals(qmmm_ff_section, "NONBONDED14%GENPOT")
      CALL section_vals_get(lj_section, n_repetition=n_lj)
      np = n_lj
      IF (n_lj /= 0) THEN
         CALL pair_potential_reallocate(inp_info%nonbonded14, 1, np, lj_charmm=.TRUE.)
         CALL read_lj_section(inp_info%nonbonded14, lj_section, start=0)
      END IF
      CALL section_vals_get(wl_section, n_repetition=n_wl)
      np = n_lj + n_wl
      IF (n_wl /= 0) THEN
         CALL pair_potential_reallocate(inp_info%nonbonded14, 1, np, williams=.TRUE.)
         CALL read_wl_section(inp_info%nonbonded14, wl_section, start=n_lj)
      END IF
      CALL section_vals_get(gd_section, n_repetition=n_gd)
      np = n_lj + n_wl + n_gd
      IF (n_gd /= 0) THEN
         CALL pair_potential_reallocate(inp_info%nonbonded14, 1, np, goodwin=.TRUE.)
         CALL read_gd_section(inp_info%nonbonded14, gd_section, start=n_lj + n_wl)
      END IF
      CALL section_vals_get(gp_section, n_repetition=n_gp)
      np = n_lj + n_wl + n_gd + n_gp
      IF (n_gp /= 0) THEN
         CALL pair_potential_reallocate(inp_info%nonbonded14, 1, np, gp=.TRUE.)
         CALL read_gp_section(inp_info%nonbonded14, gp_section, start=n_lj + n_wl + n_gd)
      END IF

   END SUBROUTINE read_qmmm_ff_section

! **************************************************************************************************
!> \brief ...
!> \param qmmm_section ...
!> \param qm_atom_index ...
!> \param qm_atom_type ...
!> \param mm_link_atoms ...
!> \param mm_link_scale_factor ...
!> \param qmmm_link ...
!> \param fist_scale_charge_link ...
!> \par History
!>      12.2004 created [tlaino]
!> \author Teodoro Laino
! **************************************************************************************************
   SUBROUTINE setup_qm_atom_list(qmmm_section, qm_atom_index, qm_atom_type, &
                                 mm_link_atoms, mm_link_scale_factor, qmmm_link, fist_scale_charge_link)
      TYPE(section_vals_type), POINTER                   :: qmmm_section
      INTEGER, DIMENSION(:), OPTIONAL, POINTER           :: qm_atom_index
      CHARACTER(len=default_string_length), &
         DIMENSION(:), OPTIONAL, POINTER                 :: qm_atom_type
      INTEGER, DIMENSION(:), OPTIONAL, POINTER           :: mm_link_atoms
      REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER     :: mm_link_scale_factor
      LOGICAL, INTENT(OUT), OPTIONAL                     :: qmmm_link
      REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER     :: fist_scale_charge_link

      CHARACTER(len=default_string_length)               :: qm_atom_kind, qm_link_element
      INTEGER                                            :: ikind, k, link_involv_mm, link_type, &
                                                            mm_index, n_var, nkind, nlinks, &
                                                            num_qm_atom_tot
      INTEGER, DIMENSION(:), POINTER                     :: mm_indexes
      LOGICAL                                            :: explicit
      REAL(KIND=dp)                                      :: scale_f
      TYPE(section_vals_type), POINTER                   :: qm_kinds, qmmm_links

      num_qm_atom_tot = 0
      link_involv_mm = 0
      nlinks = 0
      !
      ! QM_KINDS
      !
      qm_kinds => section_vals_get_subs_vals(qmmm_section, "QM_KIND")
      CALL section_vals_get(qm_kinds, n_repetition=nkind)
      DO ikind = 1, nkind
         CALL section_vals_val_get(qm_kinds, "MM_INDEX", i_rep_section=ikind, n_rep_val=n_var)
         DO k = 1, n_var
            CALL section_vals_val_get(qm_kinds, "MM_INDEX", i_rep_section=ikind, i_rep_val=k, &
                                      i_vals=mm_indexes)
            num_qm_atom_tot = num_qm_atom_tot + SIZE(mm_indexes)
         END DO
      END DO
      !
      ! QM/MM LINKS
      !
      qmmm_links => section_vals_get_subs_vals(qmmm_section, "LINK")
      CALL section_vals_get(qmmm_links, explicit=explicit)
      IF (explicit) THEN
         qmmm_link = .TRUE.
         CALL section_vals_get(qmmm_links, n_repetition=nlinks)
         ! Take care of the various link types
         DO ikind = 1, nlinks
            CALL section_vals_val_get(qmmm_links, "LINK_TYPE", i_rep_section=ikind, &
                                      i_val=link_type)
            SELECT CASE (link_type)
            CASE (do_qmmm_link_imomm)
               num_qm_atom_tot = num_qm_atom_tot + 1
               link_involv_mm = link_involv_mm + 1
            CASE (do_qmmm_link_pseudo)
               num_qm_atom_tot = num_qm_atom_tot + 1
            CASE (do_qmmm_link_gho)
               ! do nothing for the moment
            CASE DEFAULT
               CPABORT("")
            END SELECT
         END DO
      END IF
      IF (PRESENT(mm_link_scale_factor) .AND. (link_involv_mm /= 0)) &
         ALLOCATE (mm_link_scale_factor(link_involv_mm))
      IF (PRESENT(fist_scale_charge_link) .AND. (link_involv_mm /= 0)) &
         ALLOCATE (fist_scale_charge_link(link_involv_mm))
      IF (PRESENT(mm_link_atoms) .AND. (link_involv_mm /= 0)) &
         ALLOCATE (mm_link_atoms(link_involv_mm))
      IF (PRESENT(qm_atom_index)) ALLOCATE (qm_atom_index(num_qm_atom_tot))
      IF (PRESENT(qm_atom_type)) ALLOCATE (qm_atom_type(num_qm_atom_tot))
      IF (PRESENT(qm_atom_index)) qm_atom_index = 0
      IF (PRESENT(qm_atom_type)) qm_atom_type = " "
      num_qm_atom_tot = 1
      DO ikind = 1, nkind
         CALL section_vals_val_get(qm_kinds, "MM_INDEX", i_rep_section=ikind, n_rep_val=n_var)
         DO k = 1, n_var
            CALL section_vals_val_get(qm_kinds, "MM_INDEX", i_rep_section=ikind, i_rep_val=k, &
                                      i_vals=mm_indexes)
            IF (PRESENT(qm_atom_index)) THEN
               qm_atom_index(num_qm_atom_tot:num_qm_atom_tot + SIZE(mm_indexes) - 1) = mm_indexes(:)
            END IF
            IF (PRESENT(qm_atom_type)) THEN
               CALL section_vals_val_get(qm_kinds, "_SECTION_PARAMETERS_", i_rep_section=ikind, &
                                         c_val=qm_atom_kind)
               qm_atom_type(num_qm_atom_tot:num_qm_atom_tot + SIZE(mm_indexes) - 1) = qm_atom_kind
            END IF
            num_qm_atom_tot = num_qm_atom_tot + SIZE(mm_indexes)
         END DO
      END DO
      IF (PRESENT(mm_link_scale_factor) .AND. (link_involv_mm /= 0)) mm_link_scale_factor = 0.0_dp
      IF (PRESENT(fist_scale_charge_link) .AND. (link_involv_mm /= 0)) fist_scale_charge_link = 0.0_dp
      IF (PRESENT(mm_link_atoms) .AND. (link_involv_mm /= 0)) mm_link_atoms = 0
      IF (explicit) THEN
         DO ikind = 1, nlinks
            IF (PRESENT(qm_atom_type)) THEN
               CALL section_vals_val_get(qmmm_links, "QM_KIND", i_rep_section=ikind, c_val=qm_link_element)
               qm_atom_type(num_qm_atom_tot:num_qm_atom_tot) = TRIM(qm_link_element)//"_LINK"
            END IF
            IF (PRESENT(qm_atom_index)) THEN
               CALL section_vals_val_get(qmmm_links, "MM_INDEX", i_rep_section=ikind, i_val=mm_index)
               CPASSERT(ALL(qm_atom_index /= mm_index))
               qm_atom_index(num_qm_atom_tot:num_qm_atom_tot) = mm_index
               num_qm_atom_tot = num_qm_atom_tot + 1
            END IF
            IF (PRESENT(mm_link_atoms) .AND. (link_involv_mm /= 0)) THEN
               CALL section_vals_val_get(qmmm_links, "MM_INDEX", i_rep_section=ikind, i_val=mm_index)
               mm_link_atoms(ikind) = mm_index
            END IF
            IF (PRESENT(mm_link_scale_factor) .AND. (link_involv_mm /= 0)) THEN
               CALL section_vals_val_get(qmmm_links, "QMMM_SCALE_FACTOR", i_rep_section=ikind, r_val=scale_f)
               mm_link_scale_factor(ikind) = scale_f
            END IF
            IF (PRESENT(fist_scale_charge_link) .AND. (link_involv_mm /= 0)) THEN
               CALL section_vals_val_get(qmmm_links, "FIST_SCALE_FACTOR", i_rep_section=ikind, r_val=scale_f)
               fist_scale_charge_link(ikind) = scale_f
            END IF
         END DO
      END IF
      CPASSERT(num_qm_atom_tot - 1 == SIZE(qm_atom_index))

   END SUBROUTINE setup_qm_atom_list

! **************************************************************************************************
!> \brief this routine sets up all variables to treat qmmm links
!> \param qmmm_section ...
!> \param qmmm_links ...
!> \param mm_el_pot_radius ...
!> \param mm_el_pot_radius_corr ...
!> \param mm_atom_index ...
!> \param iw ...
!> \par History
!>      12.2004 created [tlaino]
!> \author Teodoro Laino
! **************************************************************************************************
   SUBROUTINE setup_qmmm_links(qmmm_section, qmmm_links, mm_el_pot_radius, mm_el_pot_radius_corr, &
                               mm_atom_index, iw)
      TYPE(section_vals_type), POINTER                   :: qmmm_section
      TYPE(qmmm_links_type), POINTER                     :: qmmm_links
      REAL(KIND=dp), DIMENSION(:), POINTER               :: mm_el_pot_radius, mm_el_pot_radius_corr
      INTEGER, DIMENSION(:), POINTER                     :: mm_atom_index
      INTEGER, INTENT(IN)                                :: iw

      INTEGER                                            :: ikind, link_type, mm_index, n_gho, &
                                                            n_imomm, n_pseudo, n_rep_val, n_tot, &
                                                            nlinks, qm_index
      INTEGER, DIMENSION(:), POINTER                     :: wrk_tmp
      REAL(KIND=dp)                                      :: alpha, my_radius
      TYPE(section_vals_type), POINTER                   :: qmmm_link_section

      NULLIFY (wrk_tmp)
      n_imomm = 0
      n_gho = 0
      n_pseudo = 0
      qmmm_link_section => section_vals_get_subs_vals(qmmm_section, "LINK")
      CALL section_vals_get(qmmm_link_section, n_repetition=nlinks)
      CPASSERT(nlinks /= 0)
      DO ikind = 1, nlinks
         CALL section_vals_val_get(qmmm_link_section, "LINK_TYPE", i_rep_section=ikind, i_val=link_type)
         IF (link_type == do_qmmm_link_imomm) n_imomm = n_imomm + 1
         IF (link_type == do_qmmm_link_gho) n_gho = n_gho + 1
         IF (link_type == do_qmmm_link_pseudo) n_pseudo = n_pseudo + 1
      END DO
      n_tot = n_imomm + n_gho + n_pseudo
      CPASSERT(n_tot /= 0)
      ALLOCATE (qmmm_links)
      NULLIFY (qmmm_links%imomm, &
               qmmm_links%pseudo)
      ! IMOMM
      IF (n_imomm /= 0) THEN
         ALLOCATE (qmmm_links%imomm(n_imomm))
         ALLOCATE (wrk_tmp(n_imomm))
         DO ikind = 1, n_imomm
            NULLIFY (qmmm_links%imomm(ikind)%link)
            ALLOCATE (qmmm_links%imomm(ikind)%link)
         END DO
         n_imomm = 0
         DO ikind = 1, nlinks
            CALL section_vals_val_get(qmmm_link_section, "LINK_TYPE", i_rep_section=ikind, i_val=link_type)
            IF (link_type == do_qmmm_link_imomm) THEN
               n_imomm = n_imomm + 1
               CALL section_vals_val_get(qmmm_link_section, "QM_INDEX", i_rep_section=ikind, i_val=qm_index)
               CALL section_vals_val_get(qmmm_link_section, "MM_INDEX", i_rep_section=ikind, i_val=mm_index)
               CALL section_vals_val_get(qmmm_link_section, "ALPHA_IMOMM", i_rep_section=ikind, r_val=alpha)
               CALL section_vals_val_get(qmmm_link_section, "RADIUS", i_rep_section=ikind, n_rep_val=n_rep_val)
               qmmm_links%imomm(n_imomm)%link%qm_index = qm_index
               qmmm_links%imomm(n_imomm)%link%mm_index = mm_index
               qmmm_links%imomm(n_imomm)%link%alpha = alpha
               wrk_tmp(n_imomm) = mm_index
               IF (n_rep_val == 1) THEN
                  CALL section_vals_val_get(qmmm_link_section, "RADIUS", i_rep_section=ikind, r_val=my_radius)
                  WHERE (mm_atom_index == mm_index) mm_el_pot_radius = my_radius
                  WHERE (mm_atom_index == mm_index) mm_el_pot_radius_corr = my_radius
               END IF
               CALL section_vals_val_get(qmmm_link_section, "CORR_RADIUS", i_rep_section=ikind, n_rep_val=n_rep_val)
               IF (n_rep_val == 1) THEN
                  CALL section_vals_val_get(qmmm_link_section, "CORR_RADIUS", i_rep_section=ikind, r_val=my_radius)
                  WHERE (mm_atom_index == mm_index) mm_el_pot_radius_corr = my_radius
               END IF
            END IF
         END DO
         !
         ! Checking the link structure
         !
         DO ikind = 1, SIZE(wrk_tmp)
            IF (COUNT(wrk_tmp == wrk_tmp(ikind)) > 1) THEN
               WRITE (iw, '(/A)') "In the IMOMM scheme no more than one QM atom can be bounded to the same MM atom."
               WRITE (iw, '(A)') "Multiple link MM atom not allowed. Check your link sections."
               CPABORT("")
            END IF
         END DO
         DEALLOCATE (wrk_tmp)
      END IF
      ! PSEUDO
      IF (n_pseudo /= 0) THEN
         ALLOCATE (qmmm_links%pseudo(n_pseudo))
         DO ikind = 1, n_pseudo
            NULLIFY (qmmm_links%pseudo(ikind)%link)
            ALLOCATE (qmmm_links%pseudo(ikind)%link)
         END DO
         n_pseudo = 0
         DO ikind = 1, nlinks
            CALL section_vals_val_get(qmmm_link_section, "LINK_TYPE", i_rep_section=ikind, i_val=link_type)
            IF (link_type == do_qmmm_link_pseudo) THEN
               n_pseudo = n_pseudo + 1
               CALL section_vals_val_get(qmmm_link_section, "QM_INDEX", i_rep_section=ikind, i_val=qm_index)
               CALL section_vals_val_get(qmmm_link_section, "MM_INDEX", i_rep_section=ikind, i_val=mm_index)
               qmmm_links%pseudo(n_pseudo)%link%qm_index = qm_index
               qmmm_links%pseudo(n_pseudo)%link%mm_index = mm_index
            END IF
         END DO
      END IF
      ! GHO
      IF (n_gho /= 0) THEN
         ! not yet implemented
         ! still to define : type, implementation into QS
         CPABORT("")
      END IF
   END SUBROUTINE setup_qmmm_links

! **************************************************************************************************
!> \brief this routine sets up all variables to treat qmmm links
!> \param qmmm_section ...
!> \param move_mm_charges ...
!> \param add_mm_charges ...
!> \param mm_atom_chrg ...
!> \param mm_el_pot_radius ...
!> \param mm_el_pot_radius_corr ...
!> \param added_charges ...
!> \param mm_atom_index ...
!> \par History
!>      12.2004 created [tlaino]
!> \author Teodoro Laino
! **************************************************************************************************
   SUBROUTINE move_or_add_atoms(qmmm_section, move_mm_charges, add_mm_charges, &
                                mm_atom_chrg, mm_el_pot_radius, mm_el_pot_radius_corr, &
                                added_charges, mm_atom_index)
      TYPE(section_vals_type), POINTER                   :: qmmm_section
      LOGICAL, INTENT(OUT)                               :: move_mm_charges, add_mm_charges
      REAL(KIND=dp), DIMENSION(:), POINTER               :: mm_atom_chrg, mm_el_pot_radius, &
                                                            mm_el_pot_radius_corr
      TYPE(add_set_type), POINTER                        :: added_charges
      INTEGER, DIMENSION(:), POINTER                     :: mm_atom_index

      INTEGER                                            :: i_add, icount, ikind, ind1, Index1, &
                                                            Index2, n_add_tot, n_adds, n_move_tot, &
                                                            n_moves, n_rep_val, nlinks
      LOGICAL                                            :: explicit
      REAL(KIND=dp)                                      :: alpha, c_radius, charge, radius
      TYPE(section_vals_type), POINTER                   :: add_section, move_section, &
                                                            qmmm_link_section

      explicit = .FALSE.
      move_mm_charges = .FALSE.
      add_mm_charges = .FALSE.
      NULLIFY (qmmm_link_section, move_section, add_section)
      qmmm_link_section => section_vals_get_subs_vals(qmmm_section, "LINK")
      CALL section_vals_get(qmmm_link_section, n_repetition=nlinks)
      CPASSERT(nlinks /= 0)
      icount = 0
      n_move_tot = 0
      n_add_tot = 0
      DO ikind = 1, nlinks
         move_section => section_vals_get_subs_vals(qmmm_link_section, "MOVE_MM_CHARGE", &
                                                    i_rep_section=ikind)
         CALL section_vals_get(move_section, n_repetition=n_moves)
         add_section => section_vals_get_subs_vals(qmmm_link_section, "ADD_MM_CHARGE", &
                                                   i_rep_section=ikind)
         CALL section_vals_get(add_section, n_repetition=n_adds)
         n_move_tot = n_move_tot + n_moves
         n_add_tot = n_add_tot + n_adds
      END DO
      icount = n_move_tot + n_add_tot
      IF (n_add_tot /= 0) add_mm_charges = .TRUE.
      IF (n_move_tot /= 0) move_mm_charges = .TRUE.
      !
      ! create add_set_type
      !
      CALL create_add_set_type(added_charges, ndim=icount)
      !
      ! Fill in structures
      !
      icount = 0
      DO ikind = 1, nlinks
         move_section => section_vals_get_subs_vals(qmmm_link_section, "MOVE_MM_CHARGE", &
                                                    i_rep_section=ikind)
         CALL section_vals_get(move_section, explicit=explicit, n_repetition=n_moves)
         !
         ! Moving charge atoms
         !
         IF (explicit) THEN
            DO i_add = 1, n_moves
               icount = icount + 1
               CALL section_vals_val_get(move_section, "ATOM_INDEX_1", i_val=Index1, i_rep_section=i_add)
               CALL section_vals_val_get(move_section, "ATOM_INDEX_2", i_val=Index2, i_rep_section=i_add)
               CALL section_vals_val_get(move_section, "ALPHA", r_val=alpha, i_rep_section=i_add)
               CALL section_vals_val_get(move_section, "RADIUS", r_val=radius, i_rep_section=i_add)
               CALL section_vals_val_get(move_section, "CORR_RADIUS", n_rep_val=n_rep_val, i_rep_section=i_add)
               c_radius = radius
               IF (n_rep_val == 1) &
                  CALL section_vals_val_get(move_section, "CORR_RADIUS", r_val=c_radius, i_rep_section=i_add)

               CALL set_add_set_type(added_charges, icount, Index1, Index2, alpha, radius, c_radius, &
                                     mm_atom_chrg=mm_atom_chrg, mm_el_pot_radius=mm_el_pot_radius, &
                                     mm_el_pot_radius_corr=mm_el_pot_radius_corr, &
                                     mm_atom_index=mm_atom_index, move=n_moves, Ind1=ind1)
            END DO
            mm_atom_chrg(ind1) = 0.0_dp
         END IF

         add_section => section_vals_get_subs_vals(qmmm_link_section, "ADD_MM_CHARGE", &
                                                   i_rep_section=ikind)
         CALL section_vals_get(add_section, explicit=explicit, n_repetition=n_adds)
         !
         ! Adding charge atoms
         !
         IF (explicit) THEN
            DO i_add = 1, n_adds
               icount = icount + 1
               CALL section_vals_val_get(add_section, "ATOM_INDEX_1", i_val=Index1, i_rep_section=i_add)
               CALL section_vals_val_get(add_section, "ATOM_INDEX_2", i_val=Index2, i_rep_section=i_add)
               CALL section_vals_val_get(add_section, "ALPHA", r_val=alpha, i_rep_section=i_add)
               CALL section_vals_val_get(add_section, "RADIUS", r_val=radius, i_rep_section=i_add)
               CALL section_vals_val_get(add_section, "CHARGE", r_val=charge, i_rep_section=i_add)
               CALL section_vals_val_get(add_section, "CORR_RADIUS", n_rep_val=n_rep_val, i_rep_section=i_add)
               c_radius = radius
               IF (n_rep_val == 1) &
                  CALL section_vals_val_get(add_section, "CORR_RADIUS", r_val=c_radius, i_rep_section=i_add)

               CALL set_add_set_type(added_charges, icount, Index1, Index2, alpha, radius, c_radius, charge, &
                                     mm_atom_chrg=mm_atom_chrg, mm_el_pot_radius=mm_el_pot_radius, &
                                     mm_el_pot_radius_corr=mm_el_pot_radius_corr, &
                                     mm_atom_index=mm_atom_index)
            END DO
         END IF
      END DO

   END SUBROUTINE move_or_add_atoms

! **************************************************************************************************
!> \brief this routine sets up all variables of the add_set_type type
!> \param added_charges ...
!> \param icount ...
!> \param Index1 ...
!> \param Index2 ...
!> \param alpha ...
!> \param radius ...
!> \param c_radius ...
!> \param charge ...
!> \param mm_atom_chrg ...
!> \param mm_el_pot_radius ...
!> \param mm_el_pot_radius_corr ...
!> \param mm_atom_index ...
!> \param move ...
!> \param ind1 ...
!> \par History
!>      12.2004 created [tlaino]
!> \author Teodoro Laino
! **************************************************************************************************
   SUBROUTINE set_add_set_type(added_charges, icount, Index1, Index2, alpha, radius, c_radius, charge, &
                               mm_atom_chrg, mm_el_pot_radius, mm_el_pot_radius_corr, mm_atom_index, move, ind1)
      TYPE(add_set_type), POINTER                        :: added_charges
      INTEGER, INTENT(IN)                                :: icount, Index1, Index2
      REAL(KIND=dp), INTENT(IN)                          :: alpha, radius, c_radius
      REAL(KIND=dp), INTENT(IN), OPTIONAL                :: charge
      REAL(KIND=dp), DIMENSION(:), POINTER               :: mm_atom_chrg, mm_el_pot_radius, &
                                                            mm_el_pot_radius_corr
      INTEGER, DIMENSION(:), POINTER                     :: mm_atom_index
      INTEGER, INTENT(in), OPTIONAL                      :: move
      INTEGER, INTENT(OUT), OPTIONAL                     :: ind1

      INTEGER                                            :: i, my_move
      REAL(KIND=dp)                                      :: my_c_radius, my_charge, my_radius

      my_move = 0
      my_radius = radius
      my_c_radius = c_radius
      IF (PRESENT(charge)) my_charge = charge
      IF (PRESENT(move)) my_move = move
      i = 1
      GetId: DO WHILE (i <= SIZE(mm_atom_index))
         IF (Index1 == mm_atom_index(i)) EXIT GetId
         i = i + 1
      END DO GetId
      IF (PRESENT(ind1)) ind1 = i
      CPASSERT(i <= SIZE(mm_atom_index))
      IF (.NOT. PRESENT(charge)) my_charge = mm_atom_chrg(i)/REAL(my_move, KIND=dp)
      IF (my_radius == 0.0_dp) my_radius = mm_el_pot_radius(i)
      IF (my_c_radius == 0.0_dp) my_c_radius = mm_el_pot_radius_corr(i)

      added_charges%add_env(icount)%Index1 = Index1
      added_charges%add_env(icount)%Index2 = Index2
      added_charges%add_env(icount)%alpha = alpha
      added_charges%mm_atom_index(icount) = icount
      added_charges%mm_atom_chrg(icount) = my_charge
      added_charges%mm_el_pot_radius(icount) = my_radius
      added_charges%mm_el_pot_radius_corr(icount) = my_c_radius
   END SUBROUTINE set_add_set_type

! **************************************************************************************************
!> \brief this routine sets up the origin of the MM cell respect to the
!>      origin of the QM cell. The origin of the QM cell is assumed to be
!>      in (0.0,0.0,0.0)...
!> \param qmmm_section ...
!> \param qmmm_env ...
!> \param qm_cell_small ...
!> \param dr ...
!> \par History
!>      02.2005 created [tlaino]
!> \author Teodoro Laino
! **************************************************************************************************
   SUBROUTINE setup_origin_mm_cell(qmmm_section, qmmm_env, qm_cell_small, &
                                   dr)
      TYPE(section_vals_type), POINTER                   :: qmmm_section
      TYPE(qmmm_env_qm_type), POINTER                    :: qmmm_env
      TYPE(cell_type), POINTER                           :: qm_cell_small
      REAL(KIND=dp), DIMENSION(3), INTENT(in)            :: dr

      LOGICAL                                            :: center_grid
      REAL(KIND=dp), DIMENSION(3)                        :: tmp
      REAL(KINd=dp), DIMENSION(:), POINTER               :: vec

! This is the vector that corrects position to apply properly the PBC

      tmp(1) = qm_cell_small%hmat(1, 1)
      tmp(2) = qm_cell_small%hmat(2, 2)
      tmp(3) = qm_cell_small%hmat(3, 3)
      CPASSERT(ALL(tmp > 0))
      qmmm_env%dOmmOqm = tmp/2.0_dp
      ! This is unit vector to translate the QM system in order to center it
      ! in QM cell
      CALL section_vals_val_get(qmmm_section, "CENTER_GRID", l_val=center_grid)
      IF (center_grid) THEN
         qmmm_env%utrasl = dr
      ELSE
         qmmm_env%utrasl = 1.0_dp
      END IF
      CALL section_vals_val_get(qmmm_section, "INITIAL_TRANSLATION_VECTOR", r_vals=vec)
      qmmm_env%transl_v = vec
   END SUBROUTINE setup_origin_mm_cell

! **************************************************************************************************
!> \brief this routine sets up list of MM atoms carrying an image charge
!> \param image_charge_section ...
!> \param qmmm_env ...
!> \param qm_atom_index ...
!> \param subsys_mm ...
!> \par History
!>      02.2012 created
!> \author Dorothea Golze
! **************************************************************************************************
   SUBROUTINE setup_image_atom_list(image_charge_section, qmmm_env, &
                                    qm_atom_index, subsys_mm)

      TYPE(section_vals_type), POINTER                   :: image_charge_section
      TYPE(qmmm_env_qm_type), POINTER                    :: qmmm_env
      INTEGER, DIMENSION(:), POINTER                     :: qm_atom_index
      TYPE(cp_subsys_type), POINTER                      :: subsys_mm

      INTEGER                                            :: atom_a, atom_b, i, j, k, max_index, &
                                                            n_var, num_const_atom, &
                                                            num_image_mm_atom
      INTEGER, DIMENSION(:), POINTER                     :: mm_indexes
      LOGICAL                                            :: fix_xyz, imageind_in_range
      TYPE(molecule_kind_type), DIMENSION(:), POINTER    :: molecule_kind

      NULLIFY (mm_indexes, molecule_kind)
      imageind_in_range = .FALSE.
      num_image_mm_atom = 0
      max_index = 0

      CALL section_vals_val_get(image_charge_section, "MM_ATOM_LIST", &
                                n_rep_val=n_var)
      DO i = 1, n_var
         CALL section_vals_val_get(image_charge_section, "MM_ATOM_LIST", &
                                   i_rep_val=i, i_vals=mm_indexes)
         num_image_mm_atom = num_image_mm_atom + SIZE(mm_indexes)
      END DO

      ALLOCATE (qmmm_env%image_charge_pot%image_mm_list(num_image_mm_atom))

      qmmm_env%image_charge_pot%image_mm_list = 0
      num_image_mm_atom = 1

      DO i = 1, n_var
         CALL section_vals_val_get(image_charge_section, "MM_ATOM_LIST", &
                                   i_rep_val=i, i_vals=mm_indexes)
         qmmm_env%image_charge_pot%image_mm_list(num_image_mm_atom:num_image_mm_atom &
                                                 + SIZE(mm_indexes) - 1) = mm_indexes(:)
         num_image_mm_atom = num_image_mm_atom + SIZE(mm_indexes)
      END DO

      ! checking, if in range, if list contains QM atoms or any atoms doubled
      num_image_mm_atom = num_image_mm_atom - 1

      max_index = SIZE(subsys_mm%particles%els)

      CPASSERT(SIZE(qmmm_env%image_charge_pot%image_mm_list) /= 0)
      imageind_in_range = (MAXVAL(qmmm_env%image_charge_pot%image_mm_list) <= max_index) &
                          .AND. (MINVAL(qmmm_env%image_charge_pot%image_mm_list) > 0)
      CPASSERT(imageind_in_range)

      DO i = 1, num_image_mm_atom
         atom_a = qmmm_env%image_charge_pot%image_mm_list(i)
         IF (ANY(qm_atom_index == atom_a)) THEN
            CPABORT("Image atom list must only contain MM atoms")
         END IF
         DO j = i + 1, num_image_mm_atom
            atom_b = qmmm_env%image_charge_pot%image_mm_list(j)
            IF (atom_a == atom_b) &
               CPABORT("There are atoms doubled in image list.")
         END DO
      END DO

      ! check if molecules in list carry constraints
      num_const_atom = 0
      fix_xyz = .TRUE.
      IF (ASSOCIATED(subsys_mm%molecule_kinds)) THEN
         IF (ASSOCIATED(subsys_mm%molecule_kinds%els)) THEN
            molecule_kind => subsys_mm%molecule_kinds%els
            DO i = 1, SIZE(molecule_kind)
               IF (.NOT. ASSOCIATED(molecule_kind(i)%fixd_list)) EXIT
               IF (.NOT. fix_xyz) EXIT
               DO j = 1, SIZE(molecule_kind(i)%fixd_list)
                  IF (.NOT. fix_xyz) EXIT
                  DO k = 1, num_image_mm_atom
                     atom_a = qmmm_env%image_charge_pot%image_mm_list(k)
                     IF (atom_a == molecule_kind(i)%fixd_list(j)%fixd) THEN
                        num_const_atom = num_const_atom + 1
                        IF (molecule_kind(i)%fixd_list(j)%itype /= use_perd_xyz) THEN
                           fix_xyz = .FALSE.
                           EXIT
                        END IF
                     END IF
                  END DO
               END DO
            END DO
         END IF
      END IF

      ! if all image atoms are constrained, calculate image matrix only
      ! once for the first MD or GEO_OPT step (for non-iterative case)
      IF (num_const_atom == num_image_mm_atom .AND. fix_xyz) THEN
         qmmm_env%image_charge_pot%state_image_matrix = calc_once
      ELSE
         qmmm_env%image_charge_pot%state_image_matrix = calc_always
      END IF

   END SUBROUTINE setup_image_atom_list

! **************************************************************************************************
!> \brief Print info on image charges
!> \param qmmm_env ...
!> \param qmmm_section ...
!> \par History
!>      03.2012 created
!> \author Dorothea Golze
! **************************************************************************************************
   SUBROUTINE print_image_charge_info(qmmm_env, qmmm_section)

      TYPE(qmmm_env_qm_type), POINTER                    :: qmmm_env
      TYPE(section_vals_type), POINTER                   :: qmmm_section

      INTEGER                                            :: iw
      REAL(KIND=dp)                                      :: eta, eta_conv, V0, V0_conv
      TYPE(cp_logger_type), POINTER                      :: logger

      logger => cp_get_default_logger()
      iw = cp_print_key_unit_nr(logger, qmmm_section, "PRINT%PROGRAM_RUN_INFO", &
                                extension=".log")
      eta = qmmm_env%image_charge_pot%eta
      eta_conv = cp_unit_from_cp2k(eta, "angstrom", power=-2)
      V0 = qmmm_env%image_charge_pot%V0
      V0_conv = cp_unit_from_cp2k(V0, "volt")

      IF (iw > 0) THEN
         WRITE (iw, FMT="(T25,A)") "IMAGE CHARGE PARAMETERS"
         WRITE (iw, FMT="(T25,A)") REPEAT("-", 23)
         WRITE (iw, FMT="(/)")
         WRITE (iw, FMT="(T2,A)") "INDEX OF MM ATOMS CARRYING AN IMAGE CHARGE:"
         WRITE (iw, FMT="(/)")

         WRITE (iw, "(7X,10I6)") qmmm_env%image_charge_pot%image_mm_list
         WRITE (iw, FMT="(/)")
         WRITE (iw, "(T2,A52,T69,F12.8)") &
            "WIDTH OF GAUSSIAN CHARGE DISTRIBUTION [angstrom^-2]:", eta_conv
         WRITE (iw, "(T2,A26,T69,F12.8)") "EXTERNAL POTENTIAL [volt]:", V0_conv
         WRITE (iw, FMT="(/,T2,A,/)") REPEAT("-", 79)
      END IF
      CALL cp_print_key_finished_output(iw, logger, qmmm_section, &
                                        "PRINT%PROGRAM_RUN_INFO")

   END SUBROUTINE print_image_charge_info

END MODULE qmmm_init
